Nontrivial stochastic resonance temperature for the kinetic Ising model 
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The kinetic Ising model in a weak oscillating magnetic field is studied in the context of stochastic 
resonance. The signal-to-noise ratio calculated with simulations is found to peak at a nontrivial res- 
onance temperature above the equilibrium critical temperature T c . We argue that its appearance is 
closely related to the vanishing of the kinetic coefficient at T c . Comparisons with various theoretical 
results in one and higher dimensions are made. 
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I. INTRODUCTION 

A series of recent papers |l|-|j| reveal unequivocally the 
phenomenon of stochastic resonance (SR) || in the ki- 
netic Ising model driven by an oscillating magnetic field. 
The possibility of SR was anticipated by viewing the Ising 
model as a system of coupled two-state oscillators in a 
stochastic force-field which is taken to be thermal fluc- 
tuations. In our previous papers we considered the 
synchronization aspects (?]] of SR manifested in the cor- 
relation function C'(T) between the magnetic field h and 
the magnetization M(T) as a function of temperature 
T. C(T) was shown to exhibit two maxima at resonance 
temperatures T r i < T c < T V 2 (T c being the critical tem- 
perature of the equilibrium system with h = 0) . The val- 
ues of both T r i and T r 2 depend on the driving frequency. 
They converge to T c in the experimentally relevant, low 
frequency limit, thus no novel characteristic temperature 
for the kinetic Ising model was observed. The results 
are qualitatively the same for two-dimensional (2D) and 
three-dimensional (3D) systems, and of course in the one- 
dimensional (ID) case T r i is absent. 

Conventionally, SR is characterized by the behavior of 
the signal-to-noise ratio (SNR). In the case of the ID 
kinetic Ising model governed by Glauber dynamics 
the SNR was computed analytically j|| both for the re- 
sponse of an embedded spin and of the whole chain, after 
Glauber's original derivation ||. Their results (to be dis- 
cussed in greater detail in Section IV) suggest that the 
SNR exhibits a maximum at a weakly frequency depen- 
dent resonance temperature T r 1D . In this manner a new 
characteristic temperature for the kinetic Ising model was 
introduced. 

The aim of this work is to study the behavior of the 
SNR for the kinetic Ising model in general spatial dimen- 
sion d. We are particularly interested in the mechanism 
responsible for the maximum of SNR versus temperature, 
its position and the associated scaling behavior as a func- 
tion of system sizes, driving amplitude, driving frequency, 
and dimensionality of the system. 



II. SIMULATION METHOD 

Standard Monte Carlo simulations for the kinetic Ising 
model in an oscillating magnetic field, h = A$m(uj s t), on 
one to four dimensional square-type lattices are carried 
out. Heat bath algorithm is used. For ID and 2D, a 
detailed study of the SNR as a function of the driving 
frequency f s = 2ttuj s , driving amplitude A, system size 
V = L d and sampling length TV are performed. Having 
determined the trend of such dependences, SNR in 3D 
and 4D are calculated only for representative values of the 
parameters. f s is expressed in units of 1/MCS (inverse 
Monte Carlo step per site) and A in units of the nearest 
neighbor coupling constant J. Free boundary conditions 
are imposed and system sizes up to V — 3 x 10 are 
simulated. 

To determine the SNR at each temperature, we fol- 
low the time evolution of the total magnetization M(t) 
for N — 2 P MCS, after sufficient equilibration, where 
p ranges from 10 to 12. The power spectrum S(f) of 
M(t) is computed using the Fast Fourier Transformation 
method. Averaging over 500 to 1000 independent runs, 
the ensemble-averaged power spectrum < S(f) > is ob- 
tained. A typical < S(f) > is presented in Fig. 1, which 
shows the characteristic sharp peak at the driving fre- 
quency f s along with background noises. To compute 
the SNR, the noise level near f s is determined by aver- 
aging < S(f) > over the interval I = [f s — 6/N,f s — 
2 / N ) \J[fs + 2/N, f s + 6/N]. The result is denoted by 
|< S(f) >|j. Taking the height of the peak minus the 
averaged noise level as the signal, we define the SNR in 
the simulations by: 



< S(f 3 ) > - \< S(f) >\i 
\<S(f)>\i 



(1) 



In a typical continuum linear-response calculation of 
the power spectrum, one obtains the following general 
form: 



S(uj) = S (uj) + Q5{uj - uj s ), 



(2) 



where Sq(oj) is the zero- field spectrum, and Q oc A is 
the amplitude of signal. These two terms correspond to 
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the noisy background and the sharp peak, respectively, as 
depicted in Fig. 1. Conventionally, the SNR is defined by 
the ratio Q/Sq(uj s ) Q. Thus, in order to compare with 
theories, some adjustments to R need to be made. It 
is straightforward to see that the proper definition after 
replacing / by lo is 

-^sim — ~J^~Rsim: 

where we have also normalized by the trivial factor N 
which arises from discrete Fourier transform so that i? S j m 
is independent of N. Hereafter our simulation results will 
be presented in terms of i? S im- Notice that in the small 
frequency and small amplitude limit, i? s im is expected to 
scale as (cf. Q) 

R sim (T;V,A) = VA 2 g(T) (4) 

where g(T) is independent of V and A. 

III. SIMULATION RESULTS 



For 3D the same trend versus frequency is observed 
with lower frequencies resulting in a higher T r , but even 
in the low-frequency limit the sharp peak is closer to T c 
than in 2D: 

T r 3D re 1.1±0.05T C , (7) 

The second peak below T c also becomes more evident for 
higher frequencies and its height is now comparable to 
that of the main peak (Fig. 5). 

Turning to 4D, namely the upper critical dimension of 
the Ising model, we see clearly a twin-peak structure in 
R(T) (Fig. 5) and the peak above T c is even closer to T c 
than in 3D: 

T, 4D w 1.05±0.05T C . (8) 

In Fig. 6 we summarize the overall trends of the height 
and position of the main resonance at T r > T c versus 
the coordination number of the lattice, 2d, in the low- 
frequency, small- amplitude limit. The decrease of the 
peak height follows roughly a power law d~ c with c«2. 



In all the dimensions considered, i? S i m (T) exhibits a 
characteristic peak at a resonance temperature T r . Re- 
sults for ID are presented in Figs. 2-3 for the effects of 
the driving frequency, system sizes and driving ampli- 
tude. From Fig. 2, we see that varying the frequency 
produces no major shift in T r . Hence we obtain 

T r 1D w J. (5) 

In Fig. 2, the scaled SNR is plotted for a variety of com- 
binations of N , V and A. The fact that they all collapse 
onto one curve confirms the expected scaling form (^) . Its 
breakdown is evident only for driving amplitudes A > J, 
as shown in Fig. 3, due to nonlinear effects. From Fig. 2 
one can also observe that for temperatures below T r the 
simulation data are sensitive to A even for quite small 
values of A. In this region much smaller A than simu- 
lated are required in order to reach the asymptotic zero- 
amplitude limit. Decreasing the amplitude further, how- 
ever, would increase the statistical errors significantly. 

The results in 2D are qualitatively the same for T > T c 
as in ID, except a stronger influence of the driving fre- 
quency on R(T) is seen (see Fig. 4). At high frequencies, 
a second peak gradually develops below T c , as illustrated 
in the magnified plot in the inset of Fig. 4. Our simula- 
tions show that this second peak becomes more refined 
as the lattice size is increased. We also check the va- 
lidity of the scaling relation (^) and find it fails only for 
fairly large driving amplitudes A > 0.2. Thus, in the low- 
frequency, small-amplitude limit, the resonance peak for 
T > T c converges to 

T 2D « 1.35±0.03T C . (6) 

With that we confirm the existence of a novel characteris- 
tic temperature distinct from T c . For higher frequencies, 
the peak above T c shifts slightly towards T c . 



IV. ANALYTICAL APPROACHES AND 
DISCUSSIONS 

A. ID exact result 

As mentioned above, the SNR for the ID kinetic Ising 
model has been computed by Schimansky-Geier et al. g 
in the low frequency limit. For completeness, their result 
is recorded here 

^ 1D = ^f^-tanh 2 (f). (9) 

Comparison with simulation data is given in Fig. 5. The 
agreement is excellent for small A, except below the peak, 
where the data is more sensitive to A and the ^4^0 
limit is not reached yet. But the overall trend of R as a 
function of A supports the above prediction. 

B. mean-field theories 

For higher dimensions, no exact result on R is available 
and we must resort to approximations. A naive approxi- 
mation is to consider at high temperature that an average 
spin in different dimensions only differ by the number of 
spins coupled to it, given by 2d. We may then make use 
of (|^), derived for ID, to obtain an approximated for- 
mula of R(T) for d > 1 by replacing J with an effective 
coupling constant 

J e = Jd. (10) 

Of course this approximation should hold only at high 
temperatures, but comparison with simulations in Fig. 5 
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(solid lines) shows agreement better than naively ex- 
pected, especially in 2D. 

Independent of the approximation, the general form of 
the SNR is given by 



fi= 7r (AM) 2 

2 S (uj s ) 



(11) 



where Sq{u)) is the frequency dependent power spectrum 
(noise strength), and AM is the amplitude of the to- 
tal magnetization induced by the external magnetic field, 
i.e., the 'signal' in M(t) = Vm + AM sin(oj s t — (f>) where 
to is the equilibrium magnetization per spin, and <f> is 
the phase shiftjf|. Note that Q = n(AM) 2 /2 is just the 
amplitude in (g). 

The simplest approach to calculate AM and Sq is the 
mean-field (MF) approximation. The mean-field AM 
can easily be found [||: 



(AM, 



MF ) — 



V 2 A 7 

J^2 



(1 -™mf) 2 " 



1 

MF 



(12) 



where tomf can be determined in the standard way by 
numerically solving the self-consistency equation tomf — 
t&nYi{2dJmyiY /T), and tmf is the mean- field relaxation 
time: 
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tmf — 



l-^(l-m^ F ) 



(13) 



The noise strength «So MF (w) can be determined using the 
Wiencr-Khintchin theorem: 



s MF M = 



2 < M 2 > MF -V 2 m 2 



MF 



tmf 



(14) 



The numerator can readily be found via the susceptibility 
to be (1 — to 2 jf )tmf- With these we obtain: 



Rmf = 



TrVA 2 (l~m 2 MF ) 
AT 2 



(15) 



We plot Rmf(T) in Fig. 5 (dashed lines). The mean-field 
result agrees with simulations at high temperatures, but 
near T c it misses the peaks entirely. Instead, the factor 
1 — mj^p yields a cusp at T™ F . Thus, mean-field theory 
fails to capture the essence of the SNR from simulations. 



C. high-temperature expansions 

The clue for the origin of the above discrepancy comes 
from re-examining the ID exact result. It is instructive 
to rewrite Rid of ([)]) in a more general form: 



R 



ID 



ttVA 2 X 
AT 



(16) 



,2J\ 



A= ^Jl-tanh 2 (^) 



T ' T cosh(2 J jT) 



(17) 



is the zero- field kinetic coefficient, given by the ratio be- 
tween the susceptibility and relaxation time 



A = 



X 



(18) 



for general dimension For ID, \ = e 2J ^ T /T and 
t = 1/[1 — tanh(2J/T)] are exact. For higher dimensions 
no exact result for x, r or A are known, but (16) remains 
valid. In fact, equation (M) is the mean-field version of 



M, with A 



MF 



1 



l MF 



)/T. 



where 



In the more refined mean-field approach introduced 
in [|| based on the time-dependent Ginzburg-Laudau 
(TDGL) equation, the SNR can be derived and it takes 
the same form as (|l6|). With that approach in 2D, 
\T is accurate up to 0(v 3 ), r up to 0(v 4 ) and hence 
AT = 1 - Av 2 + 0{v 4 ), where v = tanh( J/T) is the usual 
high-temperature expansion (HTE) parameter. This also 
yields a smooth R(T) with no peak. 

Of course the HTE of AT is available in the literature 
to much higher order for the 2D kinetic Ising model, at 
least up to 0(v 20 ) JToj] . It is interesting to ask if this 
more accurate A gives the peak in R(T). Fig. 7 plots the 
cumulative contributions up to various orders v n , and it 
is clear that even 0(v 20 ) is not enough to yield the simu- 
lated shape of R(T). That this is the right conclusion is 
indicated by doing the same procedure for ID where the 
HTE is known to all orders. The result (Fig. 7) shows the 
same trend as in 2D. The peak is not revealed until up to 
n « 70. We conclude that the peak in R(T) is an elusive 
quantity to get, its absence is due to the inaccuracy of 
our approximation of the kinetic coefficient A near T c . In 
retrospect, approximations of high-temperature nature 
are bound to fail, because their expansion parameters 
are close to unity near where the peak is supposed to be. 
We must therefore address the critical region. 



D. critical dynamics 

From the general form (|l6|), the SNR for weak field is 
proportional to the kinetic coefficient A. From renormal- 
ization group analysis ||, r ~ e~ uz , and \ ~ e _7 j where 
e oc T - T c , we get A ~ e v ( z - 2+r >), where 7 = v(2 - rj) 
is used. It is then clear that R(T) must exhibit a max- 
imum near T c because critical slowing down entails the 
vanishing of the kinetic coefficient at T c ||, R(T) must 
bend down near T c as T is lowered. 

It is, however, clear from the simulations that R(T) 
does not really vanish at T c (finite-size effect plays no role 
since the convergence with respect to V already occurs at 
rather small V) . The physical reason is that the presence 
of an oscillating magnetic field prevents the system from 
fully developing its correlations near T c within the finite 
period l// s . Singularities are rounded. In particular, 
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critical slowing down is suppressed due to the cutoff of 
t at l/f s . Hence roundings are expected to occur at 
roughly T* , where (T* - T c )/T c ~ fj vz . For infinite 
system size, the vanishing of R is then controlled by f s 
according to R(T C ) ~ fi*-2+n)l* is controlled by L 

instead if L < f s 1 ^ 2 ). Thus, for small frequencies, we 
expect the kinetic coefficient to drop significantly as T c 
is approached, giving rise to the peak in the SNR, hence 
a novel temperature scale distinct from T c . 

Since an accurate functional form of A(T) near T c is 
lacking, to partially remedy the situation we illustrate the 
above idea by means of a phcnomenological description 
for the SNR. We simply replace the above '~' signs by 
equalities and specify e = (T — T c )/T to get trg = e_l/z , 
XRGT = e-^ 2 -") and 



A 



RG 



(19) 



For the 2D kinetic Ising model, the exponent values are 
7] = 1/4, v = 1, z « 2.16 (z is not known exactly [(Tof). 
This particular form has the obvious advantage of cap- 
turing both the correct high-temperature value (both r 
and yT — > 1 as T — > oo) and the behavior near T c . In 
fact, trg and xrg agree surprisingly well with the HTE 
away from T c . The corresponding SNR, Rrg(T), is plot- 
ted in Fig. 8, which indeed shows the characteristic peak 
at about the right place. The scenario is qualitatively 
the same in higher dimensions, but the singularities are 
weaker (logarithmic in 4D) which plausibly explains why 
T r — T c decreases with d (Fig. 6). 



V. CONCLUSION 

We have simulated the kinetic Ising model in various 
spatial dimensions under the influence of an oscillating 
magnetic field. We focus on the signal-to-noise ratio, R, 
as a measure of stochastic resonance for the spins in re- 
sponse to the external field. For all the dimensions we 
study, R exhibits a clear maximum at a resonance tem- 
perature distinct from the equilibrium critical tempera- 
ture. Various theoretical approaches to calculate R are 
discussed which, when properly incorporating the critical 
slowing down at T c , agree with the simulated results. 

We have confined our attention to the case of weak 
fields, since strong fields induce more complex response 
than is treatable by linear response theories. One such 
complications is the saturation of the magnetization 
within one cycle of oscillation which would generate 
higher harmonics in the power spectrum, thus invalidat- 
ing even the usual definition of SNR. Theoretically, we 
have also confined mostly to T > T c . Below T c , localized 
excitations such as nucleation of droplets may become 
important in certain region of the parameter space. They 
are more difficult to handle than a spatially uniform per- 
turbation done here. 



With respect to other resonance temperatures defined 
by means of the correlation function between magnetiza- 
tion and external field (|]|], the present resonance tem- 
perature seems to be unrelated. Since, unlike the pre- 
vious ones, it does not converge to T c in the small fre- 
quency limit, it offers a more robust characterization of 
the phenomena of stochastic resonance in kinetic Ising 
systems. We expect that experimental measurement of 
the R(T) curve performed on mono-domain magnetic 
particles with localized magnetic moments could reveal 
this novel resonance temperature. 
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FIG. 1. Characteristic shape of < S(f) > for a driving frequency f s « 0.097. 
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FIG. 2. General shape of the R(T)/VA 2 curve in ID for different driving frequency f B , driving amplitude A, system size V, 
and time step N. The continuous line is the exact theoretical result (H). 
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FIG. 3. Break-down of the scaling law <ffl for high driving field intensity, A. JV = 4096, V = 1000 and f a = 0.019531 for all 
curves. 
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FIG. 4. General shape and frequency dependence of the R(T)/VA 2 curves in 2D. All results are for V — 50 x 50, N — 1024 
and A = 0.11. The magnified region shows the peak under T c , observable for high driving frequencies. 



G 





FIG. 6. Trend for T r and height of the peak as a function of the dimensionality of the square-type lattices. The decrease of 
the height is well fitted by a dT c power law with c = 2. 
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FIG. 8. Result using the phenomenological approximation ( |l9| ) for the kinetic coefficient (solid line) in comparison with 2D 
simulation data. 
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